WS 2016/2017
setwd("YOUR-PATH")
raw_data <- read.table(file = "data/FS_UTM.csv",
header = TRUE,
sep = ";",
stringsAsFactors = FALSE,
quote = "")
## check for issues/problems
## str(raw_data)
## apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
raw_data$settlement[raw_data$settlement > 1] <- 1
library(sp)
sites <- raw_data
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")
Please go to http://srtm.csi.cgiar.org/, download the appropriate SRTM scence(s) for our Romanian research area, unzip them and load them into R using the rgdal package.
Plot the sites onto the raster to see whether the srtm scene covers all the points.
library(rgdal)
srtm <- readGDAL("data/srtm_41_03.tif")
## data/srtm_41_03.tif has GDAL driver GTiff ## and has 6001 rows and 6001 columns
image(srtm)
points(spTransform(sites, CRS("+init=epsg:4326")))
library(raster)
srtm <- raster("data/srtm_41_03.tif")
plot(srtm)
points(spTransform(sites, CRS("+init=epsg:4326")))
## takes a bit of time
srtm <- projectRaster(srtm, res=90, crs=CRS("+init=epsg:32634"))
plot(srtm)
points(sites)
Ask the manual: where is the difference between crop and mask?
srtm <- crop(x = srtm, y = extent(sites)+20000) plot(srtm) points(sites)
Since we do not want to reproject and crop in the subsequent analyses we just export the SRTM.
writeRaster(x = srtm,
filename = "./results/srtm.tif",
overwrite = TRUE
)
Using the raster package we can create a lot of different terrain parameters by just one line of code. The result will be a multi-layer raster object, like you know it from multi-/hyperspectral satellite images. You can work on such scenes like you would work on ordinary vector objects in R.
srtm.tp <- terrain(x = srtm,
opt = c("slope",
"aspect",
"TPI",
"TRI",
"roughness",
"flowdir"),
unit = "degrees",
neighbors = 8)
plot(srtm.tp)
You should be familiar with slope and aspect but do you know the what the other parameter are? A look in the help offers insights (they are based on Wilson & Gallant (2000)):
?raster::terrain)The help also explains that we can use focal functions in order the adapt the approaches to our needs. A focal function corresponds to a moving window. It uses a matrix of weights for the neighborhood of the focal cells.
tpiw <- function(x, w=5) {
m <- matrix(1/(w^2-1), nc=w, nr=w) # mean of the surrounding pixel
m[ceiling(0.5 * length(m))] <- 0 # set the centre cell 0
f <- focal(x, m, pad = TRUE) # apply moving window
x - f
}
tpi15 <- tpiw(x = srtm, w=15)
tpi31 <- tpiw(x = srtm, w=31)
#par(mfrow=c(1,3))
#plot(srtm.tp$tpi)
#plot(tpi15)
#plot(tpi31)
For better visualisation we calculate a hillshade and overlay the TPI raster.
srtm.hs <- hillShade(slope = terrain(srtm, opt="slope"),
aspect = terrain(srtm, opt="aspect"),
angle= 150, direction = 45, normalize = TRUE)
plot(srtm.hs,
col = gray.colors(n= 255, start=.2, end=.9),
legend = FALSE)
plot(tpi31,
col=colorRampPalette(c("red", "white", "blue"))(255),
alpha = .5, add=TRUE)
writeRaster(x = tpi31,
filename = "./results/tpi31.tif",
overwrite = TRUE
)
list.files(path = "./results")
## [1] "srtm.tif" "tpi31.tif"
Now, it is time to integrate raster and points.
tmp <- extract(x = tpi31, y = sites) head(tmp)
## [1] -1.3770858 1.0698603 0.5593438 -0.6684784 2.4236159 0.8889300
length(tmp)
## [1] 508
hist(tmp)
Task: please, make it nice, i.e. visually appealing and informative.
library(rasterVis) #srtmTheme <- rasterTheme(region=terrain.colors(200)) #levelplot(srtm, par.settings = srtmTheme) levelplot(brick(tpi31, tpi15,srtm.tp$tpi), par.settings = RdBuTheme)
Wilson, J.P. & Gallant, J.C. eds., 2000. Terrain analysis: Principles and applications, New York: Wiley.